
###function to make a coefficient plot
coef_clusterplot = function(ModelResults, varnames = NULL, data, clusterid){
       library(multiwayvcov)
       library(lmtest)
       library(dplyr)
       require(ggplot2)
       library(gridExtra)
       cluster <- data[,clusterid]
       vcov_cluster <- cluster.vcov(ModelResults, cluster)
       coef_cluster <- coeftest(ModelResults, vcov = vcov_cluster)
       
       coef_cluster <- round(coef_cluster, digits = 4)
       modelcoef = coef_cluster[,"Estimate"]
       modelse = coef_cluster[,"Std. Error"]
       ylo = modelcoef - 1.96*(modelse)
       yhi = modelcoef + 1.96*(modelse)
       names = rownames(coef_cluster)
       dfplot = data.frame(names, modelcoef, modelse, ylo, yhi)
       #dfplot$names <- as.character(dfplot$names)
       dfplot <- dfplot %>%
              dplyr::mutate(color_sig = ifelse(ylo < 0 & yhi > 0, "#ef8a62", "#2c7bb6"))
       if(!is.null(varnames)){
              dfplot <- dfplot[match(varnames, dfplot$names),]
       }
       else{dfplot <- dfplot}
       
       p <- ggplot(dfplot, aes(x=names, y=modelcoef, ymin=ylo, ymax=yhi)) +
              geom_pointrange(size=1, shape=16, colour=dfplot$color_sig) +
              coord_flip() +
              geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
              xlab('') + ylab('') +
              theme(legend.position="none",
                    axis.text = element_text(size=14),
                    axis.title=element_text(size=14))
       return(p)
}

coef_comparison = function(ModelResults,n.sim = 1000, varname, data,
                           val1, val2, clusterid){
       require(arm)
       library(multiwayvcov)
       library(lmtest)
       library(gridExtra)
       library(viridis)
       library(ggridges)
       library(dplyr)
       library(ggplot2)
       library(tidyr)
       library(stringr)
       library(purrr)
       
       cluster <- data[,clusterid]
       
       vcov_cluster = lapply(ModelResults, function(x) FUN = cluster.vcov(x, cluster))
       modSumm = mapply(function(x, y) coeftest(x,y), x = ModelResults, y = vcov_cluster)
       noModels=length(modSumm)
       model_names = paste('Model',1:noModels)
       
       set.seed(12345)
       sim = Map(function(x, y) mvrnorm(n= n.sim, coef(x), y), ModelResults, vcov_cluster)
       
       fd <- list()
       
       for (i in 1:length(ModelResults)){
              X1 <- model.matrix(ModelResults[[i]])
              X2 <- model.matrix(ModelResults[[i]])
              X = model.matrix(ModelResults[[i]])
              value = as.numeric(quantile(X[, varname], c(val1, val2)))
              X1[, varname] = value[1]
              X2[, varname] = value[2]
              fd[[model_names[i]]] =  apply(apply(X2, 1, function (x) plogis(sim[[i]] %*% x)) -
                                                   apply(X1, 1, function (x) plogis(sim[[i]] %*% x)), 1, mean)
       }
       
       df_plot <- map(fd, data.frame) %>%
              map2_df(., names(.), ~mutate(.x, id = .y))
       
       names(df_plot)[1] <- "value"
       
       
       p =    ggplot(df_plot, aes(x = value, y = id, height=..density.., fill = id)) +
              geom_density_ridges(col = "grey70", scale = .8, show.legend = F) +
              scale_fill_viridis(discrete = TRUE) +
              geom_vline(xintercept = 0, colour = gray(1/2), lty = 2) +
              theme_ridges(font_size = 13, grid = F, center_axis_labels = T) +
              theme(axis.title.y = element_blank())
       return(p)
       
}

####function to make density plot based on simulation from clustering standard errors
density_clusterplot <- function(ModelResults, n.sim = 1000, varname,
                                data, clusterid, label = c("lab1", "lab2")){
       require(arm)
       library(multiwayvcov)
       library(lmtest)
       library(ggplot2)
       cluster <- data[,clusterid]
       #Cluster by ccode
       vcov_cluster <- cluster.vcov(ModelResults, cluster)
       coef_cluster <- coeftest(ModelResults, vcov = vcov_cluster)
       set.seed(12345)
       sim <- mvrnorm(n= n.sim, coef(ModelResults), vcov_cluster)
       X1 <- model.matrix(ModelResults)
       X2 <- model.matrix(ModelResults)
       X1[,varname] <- 0
       non <- apply(apply(X1, 1, function (x) plogis(sim %*% x)), 1, mean) #1 indicates row, 2= columns
       X2[,varname] <- 1
       yes <- apply(apply(X2, 1, function (x) plogis(sim %*% x)), 1, mean)
       simprodicted <- data.frame(non = non,
                                  yes = yes)
       require(reshape2)
       require(ggplot2)
       simprodicted <- melt(simprodicted)
       levels(simprodicted$variable)
       levels(simprodicted$variable) <- label
       library(ggthemes)
       legend_title <- ""
       p <- ggplot(simprodicted, aes(x=value, fill=variable)) + theme_tufte() +
              geom_density(alpha=0.4) +
              labs(x="Predicted Probability", title="",fill="") +
              theme(legend.position="bottom") + labs(x = "", caption =  paste0("Note: density plots are based on ", n.sim," MC iterations")) +
              scale_fill_manual(legend_title, values = c("gold","blue"))
       return(p)
}

###function to make a coefficient plot
density_ridgeplot = function(ModelResults, n.sim = 1000, varname, data,
                             val1, val2, clusterid){
       require(arm)
       library(multiwayvcov)
       library(lmtest)
       library(gridExtra)
       library(viridis)
       library(ggridges)
       library(dplyr)
       library(ggplot2)
       library(tidyr)
       library(stringr)
       library(purrr)
       cluster <- data[,clusterid]
       vcov_cluster <- cluster.vcov(ModelResults, cluster)
       coef_cluster <- coeftest(ModelResults, vcov = vcov_cluster)
       set.seed(12345)
       
       sim <- mvrnorm(n= n.sim, coef(ModelResults), vcov_cluster)
       
       fd <- list()
       
       for (i in 1:length(varname)){
              X1 <- model.matrix(ModelResults)
              X2 <- model.matrix(ModelResults)
              X = model.matrix(ModelResults)
              value = as.numeric(quantile(X[, varname[i]], c(val1, val2)))
              #max = as.numeric(quantile(X[, varname[i]], val2))
              X1[, varname[i]] = value[1]
              #X1[, varname[i]] = min(X1[, varname[i]])
              #X1[, varname[i]] = quantile(X1[, varname[i]], val1)
              X2[, varname[i]] = value[2]
              #X2[, varname[i]] = max(X2[, varname[i]])
              #X2[, varname[i]] = quantile(X2[, varname[i]], val2)
              
              fd[[varname[i]]] =  apply(apply(X2, 1, function (x) plogis(sim %*% x)) -
                                               apply(X1, 1, function (x) plogis(sim %*% x)), 1, mean)
       }
       
       df_plot <- map(fd, data.frame) %>%
              map2_df(., names(.), ~mutate(.x, id = .y))
       
       names(df_plot)[1] <- "value"
       
       
       p =    ggplot(df_plot, aes(x = value, y = id, height=..density.., fill = id)) +
              geom_density_ridges(col = "grey70", scale = .8, show.legend = F) +
              scale_fill_viridis(discrete = TRUE) +
              geom_vline(xintercept = 0, colour = gray(1/2), lty = 2) +
              theme_ridges(font_size = 13, grid = F, center_axis_labels = T) +
              theme(axis.title.y = element_blank())
       return(p)
}



##function for first difference
plot_roc <- function(ModelResults,linetypes = c("solid", "dotted"), interval){
       require(plotROC)
       require(pROC)
       require(dplyr)
       require(ggplot2)
       library(purrr)
       library(tibble)
       library(magrittr)
       
       pred_dv =lapply(ModelResults, function(x)
              FUN = predict(x, type = "response"))
       
       Y = lapply(ModelResults, function(x) FUN = x$y)
       
       roc = Map(function(x, y) roc(x,y), Y, pred_dv)
       
       roc_df  = lapply(roc, function(x)
              FUN= data.frame(plotx = x$specificities,
                              ploty = rev(x$sensitivities),
                              name = paste("AUC =",
                                           sprintf("%.3f",x$auc)))) %>%
              map_df(., rbind)
       
       breaks = seq(0, 1, interval)
       
       p = ggplot(roc_df, aes(x = plotx, y = ploty, color = name, linetype = name)) +
              geom_line() +
              scale_colour_discrete("") +
              scale_linetype_manual(name = '',
                                    values=linetypes) +
              geom_segment(aes(x = 0, y = 0, xend = 1, yend = 1), alpha = 0.5) +
              #geom_step(aes(linetype = name)) +
              scale_x_continuous(name = "1-Specificity",limits = c(0,1),
                                 breaks = breaks, expand = c(0.001,0.001)) +
              scale_y_continuous(name = "Sensitivity", limits = c(0,1),
                                 breaks = breaks, expand = c(0.001,0.001)) +
              theme_grey() +
              theme(legend.position=c(.8, 0.1),
                    axis.title.y = element_text(margin = margin(1,1,1,1)),
                    axis.text = element_text(size=13),
                    axis.title=element_text(size=13))
       return(p)
       
}

# Generate 1000 new datasets #
sim_cluster <- function(ModelResults, n.sim = 1000, data, clusterid){
       require(arm)
       library(multiwayvcov)
       library(lmtest)
       library(ggplot2)
       cluster <- data[,clusterid]
       vcov_cluster <- cluster.vcov(ModelResults, cluster)
       coef_cluster <- coeftest(ModelResults, vcov = vcov_cluster)
       set.seed(12345)
       sim <- mvrnorm(n= n.sim, coef(ModelResults), vcov_cluster)
       return(sim)
}



